rm(list=ls())
gc()

## ESTIMATE TSMART PID MODELS


## LOAD PACKAGES
library(data.table)
library(lfe)
library(stringr)



data = fread(paste0('panel-analysis-2008-2012.csv.gz'),select = c('DemDiff','RepDiff','DemSpExpDiff_nohh', 'RepSpExpDiff_nohh','DemSpExp_nohh_year1','RepSpExp_nohh_year1','Age_year1','Race','Gender',
                                                                           'ZipCode','countyfips','Party_year1','hh.n.adj_year1','hh.n.adj_year2',
                                                                           'hh.d.adj_year1','hh.d.adj_year2','hh.r.adj_year1','hh.r.adj_year2',
                                                                           'State', 'WhiteBlockGroupDiff','AgeBlockGroupDiff','RegsBlockGroupDiff',
                                                                           'HHIncomeBlockGroupDiff'  , 'HomeownerBlockGroupDiff' , 'YearBuiltBlockGroupDiff' , 
                                                                           'DriveWorkBlockGroupDiff'  , 'HouseValueBlockGroupDiff'))


data[!Party_year1%in%c('Democrat','Republican'), Party_year1:='Other']
data[countyfips=='',countyfips:=NA]
data[ZipCode=='',ZipCode:=NA]
data[Gender=='',Gender:=NA]
# create household change variables
data[,hh.n.diff:=hh.n.adj_year2-hh.n.adj_year1]
data[,hh.d.diff:=hh.d.adj_year2-hh.d.adj_year1]
data[,hh.r.diff:=hh.r.adj_year2-hh.r.adj_year1]


# split data
dems = data[Party_year1=='Democrat']
reps = data[Party_year1=='Republican']
oths = data[Party_year1=='Other']

rm(data)
gc()

# construct grouping variables based on age decile, gender, race, year1 treatment decile, state, year 1 household composition, and Zip Code
# we can drop files where any of these are missing

dems = dems[!is.na(Age_year1) & !is.na(Gender) & !is.na(Race)  &!is.na(ZipCode)]
reps = reps[!is.na(Age_year1) & !is.na(Gender) & !is.na(Race)& !is.na(ZipCode)]
oths = oths[!is.na(Age_year1) & !is.na(Gender) & !is.na(Race) &  !is.na(ZipCode)]

dems[,AgeDecile:=cut(Age_year1,breaks=unique(as.numeric(quantile(Age_year1, probs=seq(0,1,by=.1),na.rm=T))), include.lowest=T)]
reps[,AgeDecile:=cut(Age_year1,breaks=unique(as.numeric(quantile(Age_year1, probs=seq(0,1,by=.1),na.rm=T))), include.lowest=T)]
oths[,AgeDecile:=cut(Age_year1,breaks=unique(as.numeric(quantile(Age_year1, probs=seq(0,1,by=.1),na.rm=T))), include.lowest=T)]


# spatial seg treatment
# democrat exposure
dems[,DemSpDecile:=cut(DemSpExp_nohh_year1,breaks=unique(as.numeric(quantile(DemSpExp_nohh_year1, probs=seq(0,1,by=.1),na.rm=T))), include.lowest=T)]
reps[,DemSpDecile:=cut(DemSpExp_nohh_year1,breaks=unique(as.numeric(quantile(DemSpExp_nohh_year1, probs=seq(0,1,by=.1),na.rm=T))), include.lowest=T)]
oths[,DemSpDecile:=cut(DemSpExp_nohh_year1,breaks=unique(as.numeric(quantile(DemSpExp_nohh_year1, probs=seq(0,1,by=.1),na.rm=T))), include.lowest=T)]

# republican exposure
dems[,RepSpDecile:=cut(RepSpExp_nohh_year1,breaks=unique(as.numeric(quantile(RepSpExp_nohh_year1, probs=seq(0,1,by=.1),na.rm=T))), include.lowest=T)]
reps[,RepSpDecile:=cut(RepSpExp_nohh_year1,breaks=unique(as.numeric(quantile(RepSpExp_nohh_year1, probs=seq(0,1,by=.1),na.rm=T))), include.lowest=T)]
oths[,RepSpDecile:=cut(RepSpExp_nohh_year1,breaks=unique(as.numeric(quantile(RepSpExp_nohh_year1, probs=seq(0,1,by=.1),na.rm=T))), include.lowest=T)]



#####
# make groups
# spatial seg treatment
dems[,GroupSpDemExp:=paste(ZipCode,AgeDecile,Race, Gender, DemSpDecile,State,hh.n.adj_year1,hh.d.adj_year1, sep ='_')]
reps[,GroupSpDemExp:=paste(ZipCode,AgeDecile,Race, Gender, DemSpDecile,State,hh.n.adj_year1,hh.d.adj_year1, sep ='_')]
oths[,GroupSpDemExp:=paste(ZipCode,AgeDecile,Race, Gender, DemSpDecile,State,hh.n.adj_year1,hh.d.adj_year1, sep ='_')]

dems[,GroupSpRepExp:=paste(ZipCode,AgeDecile,Race, Gender, RepSpDecile,State,hh.n.adj_year1,hh.r.adj_year1, sep ='_')]
reps[,GroupSpRepExp:=paste(ZipCode,AgeDecile,Race, Gender, RepSpDecile,State,hh.n.adj_year1,hh.r.adj_year1, sep ='_')]
oths[,GroupSpRepExp:=paste(ZipCode,AgeDecile,Race, Gender, RepSpDecile,State,hh.n.adj_year1,hh.r.adj_year1, sep ='_')]



## output SDs


#   ###
# estimate models


dems.sd.DemSpExp = dems[!grepl('_NA_',GroupSpDemExp),list(n=.N, mean = mean(DemSpExpDiff_nohh,na.rm=T), sd = sd(DemSpExpDiff_nohh,na.rm=T)),by='GroupSpDemExp']
dems.sd.RepSpExp = dems[!grepl('_NA_',GroupSpRepExp),list(n=.N, mean = mean(RepSpExpDiff_nohh,na.rm=T), sd = sd(RepSpExpDiff_nohh,na.rm=T)),by='GroupSpRepExp']

reps.sd.DemSpExp = reps[!grepl('_NA_',GroupSpDemExp),list(n=.N, mean = mean(DemSpExpDiff_nohh,na.rm=T), sd = sd(DemSpExpDiff_nohh,na.rm=T)),by='GroupSpDemExp']
reps.sd.RepSpExp = reps[!grepl('_NA_',GroupSpRepExp),list(n=.N, mean = mean(RepSpExpDiff_nohh,na.rm=T), sd = sd(RepSpExpDiff_nohh,na.rm=T)),by='GroupSpRepExp']

oths.sd.DemSpExp = oths[!grepl('_NA_',GroupSpDemExp),list(n=.N, mean = mean(DemSpExpDiff_nohh,na.rm=T), sd = sd(DemSpExpDiff_nohh,na.rm=T)),by='GroupSpDemExp']
oths.sd.RepSpExp = oths[!grepl('_NA_',GroupSpRepExp),list(n=.N, mean = mean(RepSpExpDiff_nohh,na.rm=T), sd = sd(RepSpExpDiff_nohh,na.rm=T)),by='GroupSpRepExp']


out = list(dems.sd.DemSpExp, dems.sd.RepSpExp, reps.sd.DemSpExp,reps.sd.RepSpExp, oths.sd.DemSpExp, oths.sd.RepSpExp)
names(out) = c('Democrats - Dem Exp',
               'Democrats - Rep Exp',
               'Republicans - Dem Exp',
               'Republicans - Rep Exp',
               'Non-partisans - Dem Exp',
               'Non-partisans - Rep Exp')

save(out, file = paste0('results/sds/sds-2008-2012-main-nohh.Rdata'))
#   ###

# Estimate models
ModelDemSpExpDems = summary(felm(DemDiff ~ DemSpExpDiff_nohh + hh.d.diff+hh.n.diff+WhiteBlockGroupDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                   HHIncomeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                   DriveWorkBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpDemExp|0|countyfips, data = dems[!grepl('_NA_',GroupSpDemExp)& !is.na(countyfips)]), robust =T)
ModelDemSpExpReps = summary(felm(DemDiff ~ DemSpExpDiff_nohh+ hh.d.diff+hh.n.diff+WhiteBlockGroupDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                   HHIncomeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                   DriveWorkBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpDemExp|0|countyfips, data = reps[!grepl('_NA_',GroupSpDemExp)&!is.na(countyfips)]), robust =T)
ModelDemSpExpOths = summary(felm(DemDiff ~ DemSpExpDiff_nohh+ hh.d.diff+hh.n.diff+WhiteBlockGroupDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                   HHIncomeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                   DriveWorkBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpDemExp|0|countyfips, data = oths[!grepl('_NA_',GroupSpDemExp)&!is.na(countyfips)]), robust =T)

ModelRepSpExpDems = summary(felm(RepDiff ~ RepSpExpDiff_nohh+ hh.r.diff+hh.n.diff+WhiteBlockGroupDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                   HHIncomeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                   DriveWorkBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpRepExp|0|countyfips, data = dems[!grepl('_NA_',GroupSpRepExp)&!is.na(countyfips)]), robust =T)
ModelRepSpExpReps = summary(felm(RepDiff ~ RepSpExpDiff_nohh+ hh.r.diff+hh.n.diff+WhiteBlockGroupDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                   HHIncomeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                   DriveWorkBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpRepExp|0|countyfips, data = reps[!grepl('_NA_',GroupSpRepExp)&!is.na(countyfips)]), robust =T)
ModelRepSpExpOths = summary(felm(RepDiff ~ RepSpExpDiff_nohh+ hh.r.diff+hh.n.diff+WhiteBlockGroupDiff+AgeBlockGroupDiff+RegsBlockGroupDiff+
                                   HHIncomeBlockGroupDiff + HomeownerBlockGroupDiff + YearBuiltBlockGroupDiff + 
                                   DriveWorkBlockGroupDiff + HouseValueBlockGroupDiff|GroupSpRepExp|0|countyfips, data = oths[!grepl('_NA_',GroupSpRepExp)&!is.na(countyfips)]), robust =T)


# store output
save(ModelDemSpExpDems, ModelDemSpExpReps, ModelDemSpExpOths, 
     ModelRepSpExpDems, ModelRepSpExpReps, ModelRepSpExpOths,
     
     
     file = paste0('results/current-results-2008-2012.Rdata'))


